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This  paper  consists  of  some  notes  on  inverse  problems  of  higher 
dimensions,  in  which  the  quantities  of  interest  (local  density  and  wave 
speed)  are  functions  of  two  or  three  spatial  variables,  e.g.  p(x,z)  and 


c(x,z),  or  p(x,y,z)  and  c(x,y,z).  These  notes  are  not  designed  to  be  complete 
but  to  summarize  the  results  obtained  so  far  in  applying  layer  stripping 
techniques  to  these  problems. 


I.  Introduction 


The  subject  of  this  paper  is  the  inverse  seismic  problem  in  dimensions 
higher  than  one,  in  which  local  density  and  wave  speed  are  functions  of  more 
than  one  spatial  variable.  To  clarify  matters,  some  terminology  is  introduced. 
The  dimension  of  an  inverse  problem  is  defined  as  the  number  of  spatial 
variables  on  which  the  quantities  of  interest  (d  and  c)  depend.  Thus,  the 
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two-dimensional  (2-D)  problem  is  the  inverse  problem  of  determining  p(x,z) 


and  c(x,z)  from  surface  measurements  of  the  displacement  u(x,  z=0,  t) ,  and  the 


fhb 


three-dimensional  (3-D)  problem  is  the  inverse  problem  of  determining  |^(x,y,z  ) 
and  c(x,y,z)  from  surface  measurements  of  the  displacement  u(x,  y,  z=0,  t) . 

Note  that  the  dimension  of  a  problem  need  not  be  the  same  as  the  dimension 


of  the  medium  for  which  it  is  defined  —  a  problem  of  given  dimension  can  be 


embedded  in  a  medium  of  higher  dimension.  For  example,  the  ‘'offset  problem 


0-  p  r  <•»'»**  l\*r  I*- 

described  in.4^-is  a  1-D  problem  embedded  in  a  2-D  medium,  while  the  point- 


source  problem  of  that  same  paper  is  a  1-D  problem  embedded  in  a  3-D  medium. 

While  a  considerable  amount  of  work  has  been  done  on  the  1-D  problem, 
much  less  has  been  done  on  the  2-D  and  3-D  problems. 
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Generalizirig  1-D  results  and  techniques  to  the  2-D  and  3-D  problems  has 
proven  to  be  very  difficult,  and  applying  other  techniques  to  these  higher¬ 
dimensional  problems  has  been  contingent  on  rather  severe  assumptions.  The 
Bom  approximation  (basically  an  assumption  that  medium  parameters  vary  slowly 
with  depth)  was  used  by  Cohen  and  Bleistein  [2] ,  and  the  WKBJ  approximation 
(an  assumption  analogous  to  geometrical  optics  in  which  energy  is  assumed  to 
propagate  along  rays)  has  been  used  by  Clayton  and  Stolt  [3] .  Raz  [4]  has 
proposed  a  migration-like  technique  that  involves  a  distorted-wave  Born  model. 
Various  assumptions  are  made,  including  a  straight-ray  approximation  between 
scattering  and  observation  points.  Results  of  a  numerical  2-D  inversion  are 
presented,  and  a  3-D  procedure  proposed. 

Newton  15]  has  described  a  general  3-D  inverse  scattering  problem  solution 
that  reconstructs  a  Schrodinger  potential  from  a  scattering  amplitude  given  as 
a  function  of  energy  and  directions  of  incident  and  scattered  particles.  Solution 
of  a  generalized  Marchenko  integral  equation  is  required.  Morawetz  and 
Kriegsmann  16]  have  proposed  an  iterative  scheme  in  which  an  initial  guess  at 
a  2-D  potential  V(x,y)  is  iteratively  refined.  In  the  numerical  examples 
presented  for  a  1-D  inverse  potential  problem,  up  to  thirteen  iterations  were 
required,  and  also  some  smoothing  to  prevent  numerical  instability.  The 
computations  and  memory  required  for  2-D  inversion  are  admitted  to  be  enormous. 

The  rest  of  this  paper  can  be  divided  into  three  sections.  First,  a 
layer-stripping  algorithm  is  given  which  reconstructs  a  2-D  density  p(x,z) , 
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with  the  assumption  of  constant  wave  speed  c.  This  is  given  more  for 
illustrating  the  application  of  layer  stripping  to  higher-dimensional  problems 
than  as  a  useful  algorithm.  Next,  the  "offset  problem"  described  in  II]  is 
generalized  to  a  2-D  problem  embedded  in  a  3-D  medium,  and  a  layer-stripping 
solution  specified.  Finally,  some  thoughts  on  generalizing  1-D  results  to 
higher  dimensions  are  given,  and  some  difficulties  in  doing  this  are  discussed. 


II.  Layer-Stripping  Reconstruction  of  p(x,z) 


In  this  section  a  recursive  layer ^stripping  algorithm  is  derived  for 
solving  the  2-D  problem  with  constant  wave  speed.  In  particular,  the  density 
p(x,z)  is  reconstructed  from  surface  observation  of  pressure  p(x,  z=0,  t)  and 
medium  acceleration  w(x,  z=0,  t) .  The  wave  speed  c  is  assumed  to  be  constant 
throughout  the  medium,  and  will  in  fact  be  taken  to  be  unity  (this  amounts  to 
scaling  depth  z  by  c) .  Of  course,  this  will  not  be  a  practical  result;  however, 
it  will  illustrate  the  application  of  the  layer-stripping  idea  to  a  higher¬ 
dimensional  problem.  Assuming  a  constant  wave  speed  removes  the  problem  of 
defining  the  wavefront,  which  turns  out  to  be  the  complicating  factor  in  applying 
layer-stripping  to  higher-dimensional  problems  (see  Section  IV) . 

This  problem  was  first  formulated  and  solved  using  layer  stripping  ideas 
by  Symes  I7J.  Symes's  approach  was  to  reconstruct  the  medium  layer  by  layer  by 
solving  a  Schrodinger  equation  in  the  lateral  variable  x  to  obtain  the  lateral 
dependence  of  density  p  at  each  depth.  The  layer  stripping  solution  to  this 
problem  using  first-order  equations  which  can  easily  be  adapted  into  a  recursive 
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algorithm  is  due  to  Levy  [8] . 

The  mathematical  technique  used  to  solve  the  partial  differential  equations 
in  this  algorithm  consists  of  propagating  the  characteristic  variables  in  depth 
z,  lateral  position  x,  and  time  t.  This  technique  has  been  applied  to  the 
problem  of  the  propagation  of  axial  shear  waves  by  Achenbach  19] ,  and  to  the 
1-D  problem  by  Santosa  and  Schwetlick  [10] ;  the  form  used  here  appears  in 
Bruckstein  et  al.  [11]. 

The  acoustic  and  stress-strain  equations  for  this  problem  are 


32p/3t2  -  -  p(3w ^/Sx  +  3wz/3z) 

(1) 

3p/3x  =  -  pwx 

(2a) 

3p/3z  =*  -  pw^ 

(2b) 

where  w  and  w  are  the  horizontal  and  vertical,  respectively,  components  of 
X  z 

acceleration.  Eliminating  wx  by  substituting  (2a)  in  (1)  yields 

3w  /3z  =  (32p/3x2  -  32p/3t ?)/p  -r  (3p/3x)  (3p/3x)/p2  (3) 

z 

It  is  assumed  that  the  pressure  p  and  vertical  acceleration  w^  (hereafter 

w  will  be  replaced  by  w,  for  convenience)  have  the  forms 
z 


w(x,z,t) 

=  6(t-z)  +  wQ(x,z,t)  l(t-z) 

(4) 

p(x,z,t) 

=  Pq (x, z , t)  l(t-z) 

(5) 

where  1 ( * )  is  the  unit  step  function.  The  impulse  in  (4)  represents  the  source 
excitation;  the  step  functions  merely  represent  the  causal  natures  of  p  and  w. 

Substituting  (4)  and  (5)  in  (2b)  and  equating  the  two  impulsive  terms 
yields 

P(x,z)  -  pQ(x,z,z)  . 

Substituting  (4)  and  (5)  in  (3)  and  equating  the  two  6(*)  terms  also 
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(6) 


yields  (6) ,  so  the  equations  are  consistent.  In  addition,  equating  the 
6(*)  terms  yields  the  added  condition 

p(x,z)  =  2Op0/at)/w0lt=2.  (7) 

This  additional  condition  ensures  that  the  problem  is  not  ill-posed  and  that 
a  unique  solution  is  forthcoming. 

The  algorithm  in  differential  form  consists  of  (2b) ,  (3),  and  (6).  In 
words  it  may  be  described  as  follows: 

Assume  all  quantities  known  at  depth  z. 

Update  to  depth  z  +  A: 

(1)  Update  pressure  p  using  (2b) .  Do  for  all  x  and  t. 

(2)  Update  acceleration  w  using  (3).  Do  for  all  x  and  t. 

(3)  Obtain  the  updated  density  p  from  (6) .  Do  for  all  x. 

The  derivatives  in  (3)  could  perhaps  best  be  accomplished  by  using 
an  FFT.  The  p  and  w  updates  are  simple  replacements,  so  the  update  can  be 
done  point  by  point. 

This  algorithm  shows  how  layer-stripping  is  carried  out  in  higher- 
dimensional  problems  —  the  updates  proceed  along  the  entire  wavefront,  point 
by  point.  The  wavefront  itself  describes  the  set  of  points  currently  being 
updated.  This  example  was  chosen  to  make  the  wavefront  as  simple  as  possible 
a  flat  planar  impulse  moving  vertically  with  unit  velocity,  when  the  wave 
speed  c  is  also  varying  in  space,  the  wavefront  becomes  distorted,  and 
characterizing  it  becomes  much  more  complicated. 

Ill .  The  2-D  Offset  Problem 

In  this  section  the  1-D  offset  problem  of  [1]  is  generalized  to  two 
dimensions.  Recall  that  in  the  1-D  offset  problem  impulsive  plane  waves 


were  incident  upon  a  2-D  medium  with  1-D  material  parameter  variation,  viz. 
p(z)  and  c(z).  Running  this  experiment  twice,  at  two  different  angles  of 
incidence,  allowed  the  recovery  of  p(z)  and  c(z)  separately.  A  generalization 
of  this  experiment  will  now  allow  p(x,z)  and  c(x,z)  to  be  recovered  separately. 

The  problem  set-up  is  as  described  in  [1],  only  now  p(x,z)  and 
c(x,z)  are  functions  of  one  lateral  coordinate  as  well  as  depth,  and  the 
impulsive  plane  wave  now  has  a  normal  lying  in  the  y-z  plane,  where  y  is  the 
other  lateral  coordinate.  This  may  be  visualized  as  a  horizontal  stack  of 
identical  inhomogeneous  plates ,  with  the  normal  to  the  impulsive  plane  wave 
having  components  in  the  direction  of  the  stacking  and  in  the  direction  of 
increasing  depth. 

The  acoustic  and  stress-strain  equations  are  now 


32p/3t2  =  -  pc2 (dw^/dx  +  3Wy/3y  +  3wz/3z) 


-  pw^  =  3p/3x 


-  pwy  =  3p/3y 

-  pw  =  3p/3z 

Z 


(8) 

(9a) 

(9b) 

(90 


where  the  quantities  are  the  same  as  in  Cl)  and  (2)  and  w^  is  the  other 
acceleration  component.  Proceeding  as  in  [1],  the  fact  that  p(x,z)  and 
c(x,z)  do  not  vary  with  y  means  that  if  the  medium  is  subject  to  an  impulsive 
plane  wave  whose  Fourier  transform  for  z  <  0  (above  the  surface)  is 
e  ^  ^xX  +  +  V\  then  the  wave  number  k^  will  not  vary  with  x  and  z  either 

above  the  surface  or  below  it.  Hence  the  Fourier  transform  of  the  pressure 
takes  the  form 

p(x,y,z,w)  =  \|>(x,z,u))  e  J  yJ  =  ^(x,z,o>)  e  i  0 


(10) 


where  0^  is  the  angle  of  incidence  for  the  plane  wave  and  c^  is  the 
(homogeneous)  wave  speed  for  z  <  0  (above  the  surface) . 

Taking  Fourier  transforms  of  (8)  and  (9) ,  substituting  (10) ,  defining 

cos2  0^(x,z)  =  1  —  c(x,z)2  sin2  9^/ cQ 2  ,  (11) 

and  converting  back  to  the  time  domed n  yields,  in  perfect  analogy  to  [1] , 

(32p/3t2)  cos2  0.(x,z)  =  -  pc2  (3w  /3x  +  3w  /3z)  .  (12) 

i  x  z 

Note  that  0^(x,z)  can  be  interpreted  as  the  angle  between  the  tangent  to  the 
actual  ray  path  at  a  point  (x,y,z)  and  its  projection  on  the  x-z  plane. 

Compare  this  to  0^(z)  in  [1],  which  was  the  angle  between  the  tangent  to  the 
ray  path  at  depth  z  and  the  z-axis.  Equation  (12)  shows  that  the  problem 
has  been  reduced  from  a  2-D  problem  embedded  in  a  3-D  medium  to  a  2-D  problem 
embedded  in  a  2-D  medium. 

Since  the  partial  derivatives  in  (9a) ,  (9c) ,  and  (12)  constitute  a  gradient 
and  divergence,  respectively,  they  must  (taken  collectively)  be  independent 
of  the  choice  of  coordinates.  Thus  we  may  change  from  x  and  z  to  the  time- 
varying  curvilinear  coordinates  s  and  e,  where  s  is  normal  to  the  (2-D)  wave- 
front  and  e  is  tangent  to  it  (see  Fig.  1) .  Note  that  s  also  represents  arc 
length  along  the  projection  of  a  ray  on  the  (x,z)  plane,  and  e  represents  arc 
length  along  a  projected  wavefront.  We  have  s  =  0  along  the  surface  (z=0) , 
and  e  =  0  along  the  ray  passing  through  the  origin  (x,z)  =  (0,0).  This 
representation  will  be  important  in  the  next  section. 

Writing  (9a),  (9c),  and  (12)  in  terms  of  s  and  e  yields 

(32p/3t2)  cos2  0,(s,e)  =  -  pc2  (3w  /3s  +  3w  /3e) 
i  sc 


(13) 


(15) 


where  wg  and  wg  are  the  components  of  acceleration  in  the  appropriate 
directions.  Eliminating  wg  gives 

(32p/3t2)  cos2  0.  (s,e)  =  -  pc2(3w  /3s) +  c232p/3e2  -  (c2/p) 

X  s 

(3p/3e)  (3p/3e) 

and  defining  the  travel  times 


dT/ds  =  1/c  ts, e) 

(16) 

dx^/dx  =  cos  0i(s,e) 

,  i=lf 2 

(17) 

for  two  experiments  with  initial  angles  of  incidence  0^  and  0 ^  allows  the 
pressure  and  acceleration  to  be  written  in  the  forms 

w*(Xe,t)  =  5(t-T. )  +  w  *(X,e,t)  l(t-X.)  (18a) 

s  l  su  x 

pi(Tfe,t)  =  p*(x,e,t)  l(t-Ti)  (18b) 

where  p1  is  the  pressure  field  resulting  from  the  experiment  at  angle  of 

incidence  0 . ,  and  similarly  for  w1 . 

i  s 

Substituting  (18)  in  (15)  and  (14a)  yields 

Pg(X,e,t  =  T\)  =  pc (X,e) /cos  @i(T,e)  (19) 

which  represents  the  information  gained  from  the  first  reflection  at  T  for  all 
e.  From  (19)  (for  both  experiments)  and  (11)  c(T,e)  may  be  found,  and  then 
p(x,e)  immediately  follows. 


Equations  (14a) ,  (15) ,  (16) ,  (17) ,  and  (19)  taken  together  thus  constitute 
a  differential  algorithm  for  computing  p(T,e)  and  c(X,e),  with  the  update  taking 


place  as  an  increment  in  the  ray  path  travel  time  T.  The  algorithm  may  be 
summarized  as  follows : 


Given;  p1(T/e/t),  w1(T/e,t),  p(T,e),  c(T,e),  cos  6. (x,e) ,  T.(T,e),  i  =  1,2 

g  1  1 

Update  all  quantities  in  X. 

Each  step  is  done  pointwise  for  all  e  and  t. 


(1) 

Update  p1  s  3pX/9x  = 

-  pc  w1 
s 

(20) 

(2) 

Update  w^  •: 3w^/3x  = 

-  [ (32p^/3t2)  cos2  0_^(T,e)  -  c232p^/3e2 

+ 

(c2/p) (3p/3e)  (3pX/3e) ]/(pc) 

(21) 

(3) 

Update  -i  3x^/3x  = 

cos  G^^e) 

(22) 

(4) 

Compute  U  s  U  (X%)  = 

.2 

(p  (xte,t  =  X^J/p  (xTe,t=X^)) 

(23) 

(5) 

Compute  c  =  c(xte)  = 

2  2  1/2 

CQ  I  (U-1)/(U  sin  02  -  sin  0^  ]  ' 

(24) 

(6) 

Compute  cos  0^  :  cos 

0^ (x*  ,e)  =  [1  -  c (x+,e) 2/cq2  sin2  GJ1^2 

(25) 

(7) 

+ 

Compute  p  :  p(T  ,e) 

=  p1(T+,e,t  =  X  +)  cos  0..(x*,e)/c(x+,  e) . 

(26) 

This  algorithm  bears  a  marked  resemblance  to  the  corresponding  algorithm 
in  II] ,  and  it  is  not  difficult  to  see  why.  In  the  1-D  offset  problem  algorithm 
updates  similar  to  those  above  were  carried  out  as  the  planar  wavefront  advanced 
from  depth  z  to  depth  z  +  A.  In  the  2-D  offset  problem  the  wavefront  is  no 
longer  a  flat  plane,  but  is  described  at  time  t  by  the  equation  T(x,z)  =  t. 

Hence  the  increment  occurs  in  ray  path  travel  time  x,  which  by  definition  is  the 
same  for  all  rays,  i.e.  all  along  the  wavefront.  When  T  is  incremented,  the 


wavefront  advances  slightly,  and  information  about  the  medium  is  obtained 
from  the  first  reflection  using  (19)  ,  which  applies  all  along  the  wavefront. 

Travel  time  X  is  used  instead  of  ray  arc  length  s ,  since  the  latter  changes 
by  varying  amounts  along  the  wavefront  in  a  given  time  increment.  In  the 
1-D  problem  there  was  no  variation  of  wave  speed  c  along  the  wavefront. 

Of  course,  it  is  still  necessary  to  convert  p(T,e)  and  c(T,e)  back 
into  the  original  (x,z)  coordinates.  This  is  done  by  a  form  of  differential 
ray  tracing  or  wave  front  tracing.  Let<J<x,e)  be  the  angle  between  a  tangent 
to  the  wavefront  at  the  point  (T,e)  and  the  (horizontal)  x-axis  (see  Fig.  1). 
Clearly  the  wavefront  will  advance  locally  in  the  direction  (j)  -  90°. 

Now,  <j>  is  of  course  a  function  of  e,  unless  the  medium  is  homogeneous. 

But  <|>  changes  with  T  due  to  variation  of  the  wave  speed  c  along  the  wavefront  — 
without  such  variation,  the  wavefont  would  retain  its  shape.  This  allows  the 
derivation  of  an  update  equation  for  <$>.  From  Fig.  2,  we  have 

tan($(T  +  $T,e)  -  <f> (T ,e))  =  (c(T,e  +  6e)  -  c(x,e))5x/6e  (27) 

and  letting  fix  and  fie  go  to  zero  yields 

3<Mx,e)/3x  =  3c(x,e)/3e  .  (28) 

This  equation  is  an  update  equation  for  <j>,  since  c(x,e)  is  assumed  to  be  known 
at  X  for  all  e,  hence  3c/3e  may  be  computed  (although  this  computation  is  not  very 
stable) . 

Now,  suppose  the  coordinates  (x,z)  associated  with  the  point  (x,e)  are 
known  for  all  e.  When  X  is  incremented  by  A,  these  coordinates  will  change 


slightly,  by  amounts  ox  and  6z.  But  clearly 

Sx(.T,e)  =  c(T,e)  A  sin  (J>(T,e)  (29a) 

6z(T,e)  =  c(T,e)  A  cos  (j>CT,e)  (29b) 


This  allows  p(x,z)  and  c(x,z)  to  be  computed  recursively,  as  follows: 

Given:  c(T,e)  ,  p(T,e),  x(T,e),  z(T,e)/Cos  4»(T,e),sin  <))(T,e)  . 

Update  all  quantities  in  T.  Each  step  is  done  for  all  e. 

i 

(1)  Update  cos  <P  from  B  cos  4>/3t  =  -  (sin  <J>)  Bc/Be  (30) 

(2)  Update  sin  <j>  from  B  sin  4>/3t  =  (cos  <)>)  3c/3e  (31) 

(3)  Update  x  and  z  from  (29) 

(4)  Update  c(T,e)  and  P(T,e)  by  the  algorithm  (18)  -  (24) 

(5)  Output  c(T  +  A,e) ,  p(T  +  A,e),  x(T  +  A,e),z(T  +  A,e)  as  c(x,z) 
and  p(x,z) .  This  is  quite  suitable  for  plotting. 

Note  that  (28)  has  been  used  in  (30)  and  (31)  ,  and  that  <P(0,e)  is  initialized 

to  zero. 

IV.  Generalizations  of  1-D  Results  to  Higher  Dimensions 

In  this  section  some  ways  in  which  methods  and  results  for  the  1-D 
inverse  seismic  problem  generalize  to  higher  dimensions  are  discussed.  Also, 
some  difficulties  in  applying  layer-stripping  to  higher-dimensional  problems 
are  noted. 

It  is  known,  e.g.  [10J ,  that  for  the  1-D  inverse  seismic  problem  in  which 
an  impulsive  plane  wave  is  normally  incident  on  a  1-D  medium,  and  the  upgoing 
wave  at  the  surface  measured,  then  the  only  information  about  the  medium  that 


can  be  reconstructed  exactly  is  the  impedance  as  a  function  of  travel 
time,  viz.  pc(i).  How  might  this  result  generalize  to  higher  dimensions? 

Writing  the  2-D  acoustic  and  stress-strain  equations  in  the  (s,e) 
coordinates  defined  in  the  last  section  gives 
2 

p  =  -  pc  (3u  /3s  +  3u  /3e)  (32) 

s  e 

3p/3s  =  -  p  32  ug/3t2  (33a) 

3p/3e  =  -  p  32  u  /3t2  (33b) 

e 

where  ug  and  u^  are  components  of  displacement  in  the  appropriate  directions. 
Now,  in  the  1-D  case  changing  variables  from  depth  to  travel  time  resulted 
in  a  set  of  equations  entirely  in  terms  of  the  impedance  pc(T),  which  allowed 
recovery  of  this  quantity  by  layer-stripping.  Unfortunately,  this  is  not 
possible  for  (32)  and  (33) ,  since  e  would  also  have  to  be  differentially  scaled 
by  c,  and  this  brings  in  other  terms.  And  as  long  as  p  and  c  are  present 
separately  in  these  equations,  there  is  no  way  they  can  be  propagated  from 
knowledge  (from  the  first  reflection)  of  their  product  Pc  alone. 

The  solution  here  is  to  recognize  an  implicit  feature  of  the  1-D  inverse 
seismic  problem:  Since  the  problem  takes  place  along  a  single  vertical  ray 
path,  only  acoustic  (i.e.  P-wave)  wave  propagation  along  this  path  need  be 
considered.  In  the  2-D  case,  this  is  tantamount  to  considering  only  acoustic 
wave  propagation  along  a  ray  path.  From  the  nature  of  acoustic  wave  propagation 
this  means  that  ug  is  negligible  ([12],  p.  95).  (Note  that  this  assumption 
would  simplify  the  algorithms  of  the  preceding  sections.)  With  this 
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assumption,  (32)  and  (33)  become 


p  =  -  pc  3u  /3s 


3p/3s  ■  -  p32ug/3t2 


(34) 

(35) 


which  have  the  same  form  as  the  basic  1-D  equations.  Defining  outgoing  and 
incoming  waves  as 

0(s,e,t)  =  p//pc  +  /pc  3u  /3t 


I(s,e,t)  =  p//pc  ~  /pc  3u  /3t 


(36a) 

(36b) 


and  assuming  an  impulse  present  in  the  outgoing  wave  yields,  as  in  [1],  the 
fast  Cholesky  equations  of  the  1-D  problem 


(3/3t  +  3/3t)  0(T,e,t)  =  -  r(x,e)  l(x,e,t) 

(3/3x  -  3/3t)  I (x,e,t)  «  -  r(x,e)  0(x,e,t) 
r  (X  ,e)  =  21  (x  ,e,X) 


(37a) 

(37b) 

(38) 


now  applied  along  each  ray  (i.e.  for  each  e) .  Thus  instead  of  reconstructing 
pc(X),  we  now  reconstruct  pc(x,e). 

A  variation  on  the  1-D  problem  provides  for  pure  shear  wave  propagation, 

with  pc (X )  again  being  reconstructed.  For  the  2-D  problem,  we  simply  neglect 

u  instead  of  u  .  Since  (32)  and  (33)  are  symmetric  in  u  and  u  ,  the  result 
s  e  s  e 

is  once  again  a  fast  Cholesky  algorithm  which  reconstructs  pc(x,e). 

As  in  the  1-D  problem,  some  sort  of  offset  experiment,  involving  the  medium 
responses  to  impulsive  plane  waves  at  two  different  angles  of  incidence,  is 
necessary  in  order  to  reconstruct  p  and  c  separately,  and  as  functions  of  x  and 
z.  The  2-D  offset  problem  where  the  normal  to  the  plane  wave  lies  in  the 
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(y,z)  plane  was  solved  in  the  previous  section.  More  desirable  would  be  a 
solution  to  the  2-D  offset  problem  where  the  normal  to  the  plane  wave  lies  in 
the  (x,z)  plane  (so  that  all  of  the  action  takes  place  in  this  plane) ,  but 
there  seems  to  be  no  way  to  relate  the  different  wave  front  histories  resulting 
from  the  two  experiments  to  each  other. 

Why  is  a  2-D  offset  experiment  necessary  in  order  to  reconstruct  p(x,z) 
and  c(x,z)  separately?  Considering  the  numbers  of  observed  quantities,  desired 
quantities,  and  dependent  variables  of  both  sheds  some  light  on  this.  The 

following  table  summarizes  the  situation  for  each  problem: 

Total  number  Total  number 

of  dependent  of  dependent 


Problem 

Observations 

variables 

i  Output  variables 

1-D 

upgoing  wave  U(0,t) 

IV 

pc  (t) 

1 

1-D  offset 

1  2 

velocities  v  (0,t),  v  (0,t) 

2 

p(z),  c(z) 

2 

2-D 

upgoing  wave  U(x,0,t) 

2 

pc (T,e) 

2 

2-D  offset 

1  2 

velocities  v  (x,0,t),  v  (x,0,t) 

4 

p (x , z) ,  c (x , z) 

4 

3-D 

upgoing  wave  U(x,y,0,t) 

3 

pc(x,e1,e2) 

3 

1-D  elastic 

upgoing  P  and  S  waves  from 

P  and  S  excitations 

3 

p(z) ,  X(z) ,  y(z) 

3 

p(x,y,z) 

acceleration  w(x,y,0,t) 

3 

p(x,y,z) 

3 

For  the  elastic  problem,  the  upgoing  waves  are  constructed  from  different 
velocity  components,  and  the  converted  reflection  responses  (P  to  S)  and 
S  to  P)  are  the  same.  Hence  there  are  only  three  measured  quantities,  instead 


of  four. 
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